# -*- coding: utf-8 -*-
"""
For NDVI Drought (1.1.24)
    

Created on Sat Jul  4 07:23:28 2020

@author: Farman
"""
import sys
import pathlib
 
base_dir = pathlib.Path(sys.argv[0]).absolute().parent.parent

if sys.path.count(base_dir) == 0:
    sys.path.append(str(base_dir))
    
    
from Auxes       import MakeDeepDir, ListFiles
from CommonBase  import ModisHdfFileMan
from CommonBase  import HdfMan
import numpy
import pyhdf.SD
import operator



'''
About VI_Quality layer of MOD13Q1

    Bit          Description
    0-1      NDVI quality
             00 VI produced, good quality;
             01 VI produced, but check other QA
             10 (pixel not produced due to cloud effects)
    				 		produced but most likley cloudy pixel
             11 pixel not produced due to other reasons than clouds
    2-5      VI usefulness
             four bit range 0= highest quality 13 quality so low that it is not useful, 
    	    14 L1B data faulty, 15 not useful for any other reason/not processed
    6-7      aerosol quantity:
              climatology (00), low (01),  average (10), high (11) 
    8        Adjacent cloud detected; (1)  YES,   (0)  NO
    9        YES (1) NO (0) atmosphere BRDF correction performed
    10       mixed clouds YES (1) NO (0)mixed clouds
    11-13    Land/Water flag  
    		000 ocean; 
    		001: Land (Nothing else but land) 
    		010: Ocean coastlines and lake shorelines 
    		011: Shallow inland water 
    		100: Ephemeral water 
    		101: Deep inland water 
    		110: Moderate or continental ocean 
    		111: Deep ocean 
    14       0 NO; 1 YES possible snow/ice
    15       0 NO; 1 YES possible shadow
    
    if filter_by_vi_quality:
        ndvi_quality, ndvi_quality_mask = get_dataframe(hdf_obj, '250m 16 days VI Quality')
        ndvi_quality = (ndvi_quality & 0x03) == 0 # NDVI quality
        ndvi_mask *= ndvi_quality
    #print('%8.5f %8.5f'%(sum(sum(ndvi_quality))/4800/48, sum(sum(ndvi_mask))/4800/48)) # percentage of good quality pixels
    # other quality check insert here.'''


#-----------------------------------------------------------------------------
# all hdf_obj is got by hdf_obj = pyhdf.SD.SD(hdf_file)

def get_dataframe(hdf_obj, frame_name):
    dataset = hdf_obj.select(frame_name)
    data_frame = dataset.get()
    fillvalue = dataset.getfillvalue()
    data_frame_mask = data_frame != fillvalue
    return data_frame, data_frame_mask


def get_dataframe_calibrated(hdf_obj, frame_name):
    '''
    hdf.select('').getcalc() :
    5-element tuple holding:
        - cal: calibration factor (attribute 'scale_factor')
        - cal_error : calibration factor error
                      (attribute 'scale_factor_err')
        - offset: calibration offset (attribute 'add_offset')
        - offset_err : offset error (attribute 'add_offset_err')
        - data_type : type of the data resulting from applying
                      the calibration formula to the dataset values
                      (attribute 'calibrated_nt')
 
    The calibration coefficients are part of the so-called
        "standard" SDS attributes. The values inside the tuple returned
        by 'getcal' are those of the following attributes, in order::
       
        scale_factor, scale_factor_err, add_offset, add_offset_err,
        calibrated_nt
 
    Original dataset values 'orival' are converted to calibrated
    values 'calval' through the formula:       
           calval = cal * (orival - offset)           
    '''
    dataset = hdf_obj.select(frame_name)
    data_frame = dataset.get()
    fillvalue = dataset.getfillvalue()
    data_frame_mask = data_frame != fillvalue
    scale_factor, scale_factor_err, add_offset, add_offset_err, calibrated_nt = dataset.getcal()
    #print(scale_factor, scale_factor_err, add_offset, add_offset_err, calibrated_nt)
    data_frame = data_frame / scale_factor + add_offset
    return data_frame * data_frame_mask, data_frame_mask


def get_dataframe_zero_as_invalid(hdf_obj, frame_name):
    dataset = hdf_obj.select(frame_name)
    data_frame = dataset.get()
    fillvalue = dataset.getfillvalue()
    data_frame_mask = data_frame != fillvalue
    data_frame_non_zero = data_frame != 0
    data_frame_mask *= data_frame_non_zero
    data_frame *= data_frame_mask
    return data_frame, data_frame_mask


def get_dataframe_calibrated_zero_as_invalid(hdf_obj, frame_name):
    dataset = hdf_obj.select(frame_name)
    data_frame = dataset.get()
    fillvalue = dataset.getfillvalue()
    data_frame_mask = data_frame != fillvalue
    data_frame_non_zero = data_frame != 0
    data_frame_mask *= data_frame_non_zero
    scale_factor, scale_factor_err, add_offset, add_offset_err, calibrated_nt = dataset.getcal()
    data_frame = data_frame / scale_factor + add_offset    
    data_frame *= data_frame_mask
    return data_frame, data_frame_mask


#-----------------------------------------------------------------------------
# For statictics

def update_min_max(val_min, val_max, val_stat_mask, val, val_mask):
    '''
    '''
    val_hot    = val_mask        # val   valid mask
    val_cool   = operator.xor(True, val_mask) # ndvi invalid mask
    stat_hot   = val_stat_mask   # stat   valid mask
    stat_cool  = operator.xor(True, val_stat_mask) # stat invalid mask
    
    # both hot region
    both_hot = val_hot * stat_hot # ndvi and stat both valid        
    both_hot_min  = (val * (val < val_min) + val_min * (val_min <= val)) * both_hot
    both_hot_max  = (val * (val > val_max) + val_max * (val_max >= val)) * both_hot

    # input is hot, stat is cool
    part_hot = val_hot * stat_cool
    part_hot_min = val * part_hot
    part_hot_max = val * part_hot
    
    #           input ndvi invalid   |      input ndvi is valid      |  input ndvi and stat are
    #                                |   but ndvi stat is invalid    |        both valid
    #-------------------------------------------------------------------------------------------------
    #                                |           part hot            |         both hot
    val_min  = val_min * val_cool    +         part_hot_min          +       both_hot_min
    val_max  = val_max * val_cool    +         part_hot_max          +       both_hot_max
    
    val_stat_mask = val_mask | val_stat_mask
    return val_min, val_max, val_stat_mask


def update_sum_count(val_sum, val_mask_count, val, val_mask):
    val_sum        += val * val_mask
    val_mask_count += 1 * val_mask
    return val_sum, val_mask_count


#-----------------------------------------------------------------------------
def calc_ndvi(nir, nir_mask, red, red_mask):
    '''
    NDVI - Normalized Difference Vegetation Index
    
           NIR - RED
    NDVI = ---------
           NIR + RED
    '''
    ndvi_mask = nir_mask * red_mask
    
    sub = nir - red
    add = nir + red
    
    mask_add = add == 0.0 # avoid divide by zero errors.
    
    ndvi = sub / (add + mask_add)
    ndvi = ndvi * ndvi_mask
    
    return ndvi, ndvi_mask # now invalid value is zero.



def get_ndvi_MOD13Q1(hdf_obj):
    '''
    Get data_frame of name frame_name from modis file given.

    Parameters
    ----------
    hdf_file : string
        name of hdf file.
    
    frame_name : string
        name of data frame in modis hdf file.
        
    Return
    ------
        numpy array. data frame get.
    '''
    nir, nir_mask = get_dataframe(hdf_obj, '250m 16 days NIR reflectance')
    red, red_mask = get_dataframe(hdf_obj, '250m 16 days red reflectance')    
    ndvi, ndvi_mask = calc_ndvi(nir, nir_mask, red, red_mask)
    return ndvi, ndvi_mask




def get_ndvi_stat(modis_files):
    ndvi_min, ndvi_max, ndvi_stat_mask = None, None, None
    ndvi_sum, ndvi_mask_count = None, None    
    
    for n in range(len(modis_files)):
        file = modis_files[n]
        
        print('Stat ndvi %2d / %2d : %s'%(n+1, len(modis_files), file))
        
        hdf_obj = pyhdf.SD.SD(file)        
        ndvi, ndvi_mask = get_ndvi_MOD13Q1(hdf_obj)
        
        if ndvi_min is None:
            ndvi_min = ndvi # numpy.resize(numpy.array( 2), ndvi.shape)
            ndvi_max = ndvi #numpy.resize(numpy.array(-2), ndvi.shape)
            ndvi_stat_mask = ndvi_mask
            ndvi_sum = ndvi
            ndvi_mask_count = 1 * ndvi_mask
        else:
            ndvi_min, ndvi_max, ndvi_stat_mask = update_min_max(ndvi_min, ndvi_max, ndvi_stat_mask, ndvi, ndvi_mask)
            ndvi_sum, ndvi_mask_count = update_sum_count(ndvi_sum, ndvi_mask_count, ndvi, ndvi_mask)
        
        hdf_obj.end()
    
    ndvi_invalid_sum_mask = ndvi_mask_count == 0
    ndvi_mean = ndvi_sum / (ndvi_mask_count + ndvi_invalid_sum_mask)
    ndvi_mean = ndvi_mean * (ndvi_mask_count != 0)
    return ndvi_min, ndvi_mean, ndvi_max, ndvi_stat_mask, ndvi_mask_count
    

#-----------------------------------------------------------------------------
     
def calc_vci(ndvi, ndvi_mask, ndvi_min, ndvi_max, ndvi_stat_mask):
    '''
    Get VCI of modis hdf file given.

    VCI - Vegetation Condition Index
    
            ndvi     - ndvi_min
    vci = ----------------------
            ndvi_max - ndvi_min
    '''
    vci_mask = ndvi_mask * ndvi_stat_mask
    
    up   = ndvi     - ndvi_min
    down = ndvi_max - ndvi_min
    
    fix_down = down == 0.0 # avoid divide by zero errors.
    mask_down =  down != 0.0
    vci_mask *= mask_down
    
    vci = up / (down + fix_down)
    vci = vci * vci_mask
    return vci, vci_mask



#-----------------------------------------------------------------------------


def calc_bt(radiation, radiation_mask, k1, k2):
    '''
    Calculate brightness temperature from radiation.
    
                  k2
    bt = ---------------------
                   k1
            ln(---------- + 1)
                L_lamda
         
    k1  k2 - constants,定标常数
    L_lamda - 热红外波段的辐亮度
    
    Recommand for test:
    k1 = 60.776
    k2 = 1560.56
    '''
    invalid = radiation_mask == False
    #return invalid, invalid
    radiation *= radiation_mask  
    radiation += invalid    
    bt = k2 / numpy.log(k1/radiation + 1) * radiation_mask
    bt_mask = radiation_mask
    return bt * bt_mask, bt_mask


def get_bt_stat(modis_files, k1, k2):
    bt_min, bt_max, bt_stat_mask = None, None, None
    bt_sum, bt_mask_count = None, None
    
    for n in range(len(modis_files)):
        file = modis_files[n]
        
        print('Stat   bt %2d / %2d : %s'%(n+1, len(modis_files), file))
        
        hdf_obj = pyhdf.SD.SD(file)     
        mir, mir_mask = get_dataframe_calibrated_zero_as_invalid(hdf_obj, '250m 16 days MIR reflectance')
        bt, bt_mask = calc_bt(mir, mir_mask, k1, k2)
        
        if bt_min is None:
            bt_min = bt # numpy.resize(numpy.array( 2), ndvi.shape)
            bt_max = bt # numpy.resize(numpy.array(-2), ndvi.shape)
            bt_stat_mask = bt_mask
            bt_sum = bt
            bt_mask_count = 1 * bt_mask
        else:
            bt_min, bt_max, bt_stat_mask = update_min_max(bt_min, bt_max, bt_stat_mask, bt, bt_mask)
            bt_sum, bt_mask_count = update_sum_count(bt_sum, bt_mask_count, bt, bt_mask)
        
        hdf_obj.end()
    
    bt_invalid_sum_mask = bt_mask_count == 0
    bt_mean = bt_sum / (bt_mask_count + bt_invalid_sum_mask)
    bt_mean = bt_mean * (bt_mask_count != 0)
    return bt_min, bt_mean, bt_max, bt_stat_mask, bt_mask_count

#-----------------------------------------------------------------------------

def calc_tci(bt, bt_mask, bt_min, bt_max, bt_stat_mask):
    '''
    Calculate Temperature Condition Index, TCI,温度状况指数
    
          BT_max - BT
    TCI = ---------------
          BT_max - BT_min
          
    Input:
        bt_list - brightness temperature data frames,
                  each frame is in type of numpy array
    Output:
        temperature condition index.
    '''
    tci_mask = bt_mask * bt_stat_mask
    
    up   = bt_max - bt
    down = bt_max - bt_min
    
    fix_down = down == 0.0 # avoid divide by zero errors.
    mask_down =  down != 0.0
    tci_mask *= mask_down
    
    tci = up / (down + fix_down)
    tci = tci * tci_mask
    return tci, tci_mask



def calc_vhi(vci, vci_mask, k_vci, tci, tci_mask, k_tci):
    '''
    Calculate vegetation healthy index, VHI, 植被健康度指数
    
    VHI = vci * k_vci + tci * k_tci
    
    Input:
        vci   - VCI data frame
        k_vci - coef. of vci
        tci   - TCI data frame
        k_tci - coef. of tci
    Output:
        vegetation healthy index.
    '''
    vhi = vci * k_vci + tci * k_tci
    vhi_mask = vci_mask * tci_mask
    return vhi, vhi_mask



def get_vhi(modis_file, 
            k1, k2, k_tci, k_vci, 
            ndvi_min, ndvi_max, ndvi_stat_mask, 
            bt_min, bt_max, bt_stat_mask):
    hdf_obj = pyhdf.SD.SD(modis_file)
    ndvi, ndvi_mask = get_ndvi_MOD13Q1(hdf_obj)        
    vci, vci_mask = calc_vci(ndvi, ndvi_mask, ndvi_min, ndvi_max, ndvi_stat_mask)
    
    mir, mir_list = get_dataframe_calibrated_zero_as_invalid(hdf_obj, '250m 16 days MIR reflectance')
    bt, bt_mask = calc_bt(mir, mir_list, k1, k2)
    tci, tci_mask = calc_tci(bt, bt_mask, bt_min, bt_max, bt_stat_mask)
    
    vhi, vhi_mask = calc_vhi(vci, vci_mask, k_vci, tci, tci_mask, k_tci)
    hdf_obj.end()
    return vhi, vhi_mask



def get_vhi_stat(modis_files, 
                 k1, k2, k_tci, k_vci, 
                 ndvi_min, ndvi_max, ndvi_stat_mask, 
                 bt_min, bt_max, bt_stat_mask):
    
    vhi_min, vhi_max, vhi_stat_mask = None, None, None
    vhi_sum, vhi_mask_count = None, None    
    
    for n in range(len(modis_files)):
        file = modis_files[n]
        print('Stat  vhi %2d / %2d : %s'%(n+1, len(modis_files), file))
        hdf_obj = pyhdf.SD.SD(file)
        ndvi, ndvi_mask = get_ndvi_MOD13Q1(hdf_obj)        
        vci, vci_mask = calc_vci(ndvi, ndvi_mask, ndvi_min, ndvi_max, ndvi_stat_mask)
        
        mir, mir_list = get_dataframe_calibrated(hdf_obj, '250m 16 days MIR reflectance')
        bt, bt_mask = calc_bt(mir, mir_list, k1, k2)
        tci, tci_mask = calc_tci(bt, bt_mask, bt_min, bt_max, bt_stat_mask)
        
        vhi, vhi_mask = calc_vhi(vci, vci_mask, k_vci, tci, tci_mask, k_tci)

        if vhi_min is None:
            vhi_min = vhi # numpy.resize(numpy.array( 2), ndvi.shape)
            vhi_max = vhi # numpy.resize(numpy.array(-2), ndvi.shape)
            vhi_stat_mask = vhi_mask
            vhi_sum = vhi
            vhi_mask_count = 1 * vhi_mask
        else:
            vhi_min, vhi_max, vhi_stat_mask = update_min_max(vhi_min, vhi_max, vhi_stat_mask, vhi, vhi_mask)
            vhi_sum, vhi_mask_count = update_sum_count(vhi_sum, vhi_mask_count, vhi, vhi_mask)
        
        hdf_obj.end()
    
    vhi_invalid_sum_mask = vhi_mask_count == 0
    vhi_mean = vhi_sum / (vhi_mask_count + vhi_invalid_sum_mask)
    vhi_mean = vhi_mean * (vhi_mask_count != 0)
    return vhi_min, vhi_mean, vhi_max, vhi_stat_mask, vhi_mask_count


def calc_vhi_drought(vhi, vhi_mask, vhi_mean, vhi_stat_mask):
    '''
    See calc_ndvi_drought.
    '''
    vhi_centered = vhi - vhi_mean
    vhi_centered_mask = vhi_mask * vhi_stat_mask
    return vhi_centered, vhi_centered_mask


def calc_vhi_drought_grade(vhi, vhi_mask, grades, vhi_min, vhi_mean, vhi_stat_mask):
    '''
    '''
    vhi_low_range = vhi_min - vhi_mean
    vhi_seg_range = vhi_low_range / (grades - 1)
    
    class_mask  = vhi_mask * vhi_stat_mask
    not_both_invalid = True ^ class_mask
    
    vhi_centered = vhi - vhi_mean
    
    classx = vhi_centered / (vhi_seg_range * class_mask + not_both_invalid)
    classx *= (classx > 0.0) * class_mask
    classx = numpy.array(classx, dtype='int')
    return classx, class_mask
    
    

#-----------------------------------------------------------------------------

'''
modis_path = '/data/yule/disk2/MOD13Q1'
modis_hdf_file_list_file = 'file_list.txt'
#ListFiles.list_files_recursive_to_file(modis_path, '*.hdf', modis_hdf_file_list_file)

# init
print('Init input varables')
modis_hdf_file_list_file = 'file_list.txt'
product = 'MOD13Q1'
year = 2015
day_of_year = 177
hv_pos = 'h27v05'
filter_by_vi_quality = False
k1 = 60.776
k2 = 1560.56
k_vci = 1
k_tci = 1
grades = 5

# get files to process
print('Get MOD13Q1 files to process')
mm = ModisHdfFileMan.ModisMan(modis_hdf_file_list_file)
modis_files = mm.select(product, '*', day_of_year, hv_pos)
modis_file = mm.select(product, year, day_of_year, hv_pos)[0]
ndvi_min, ndvi_mean, ndvi_max, ndvi_stat_mask, ndvi_mask_count = get_ndvi_stat(modis_files)
bt_min, bt_mean, bt_max, bt_stat_mask, bt_mask_count = get_bt_stat(modis_files, k1, k2)
vhi_min, vhi_mean, vhi_max, vhi_stat_mask, vhi_mask_count = get_vhi_stat(modis_files, 
                                                                         k1, k2, k_tci, k_vci, 
                                                                         ndvi_min, ndvi_max, ndvi_stat_mask, 
                                                                         bt_min, bt_max, bt_stat_mask)
vhi, vhi_mask = get_vhi(modis_file, k1, k2, k_tci, k_vci, ndvi_min, ndvi_max, ndvi_stat_mask, bt_min, bt_max, bt_stat_mask)


'''




#-----------------------------------------------------------------------------

def calculate_vhi_drought_with_coor(modis_hdf_file_list_file, 
                                 hv_pos, 
                                 year, 
                                 day_of_year,
                                 k1, k2, k_vci, k_tci,
                                 grades,                                 
                                 drought_save_path,
                                 filter_by_vi_quality):
    '''
    Get vhi drought value of position hv_pos.
    
    hv_pos : hv_position or hv_position list.
    year and day_of_year should NOT be "*".
    
    vhi drought value is saved to '250m 16 days NDVI' layer of dest hdf file.
    vhi drought grade is saved to '250m 16 days VI Quality' layer of dest hdf file.
    
    ndvi drought value saved is 10000 times of actual.
    
    filter_by_vi_quality : ignored now. see code.
    '''
    filter_by_vi_quality = int(filter_by_vi_quality)
    year, day_of_year, grades = int(year), int(day_of_year), int(grades)
    k1, k2, k_vci, k_tci = float(k1), float(k2), float(k_vci), float(k_tci)
    product = 'MOD13Q1'
    
    # prepare path for result saving.
    MakeDeepDir.make_path(drought_save_path)

    # prepare modis file manager
    mm = ModisHdfFileMan.ModisMan(modis_hdf_file_list_file)
    
    # convert string hv_pos to list, for elastic process of multi positions.
    if type(hv_pos) == str:
        hv_pos = [hv_pos]
    
    hv_positions = hv_pos
    
    # process for each hv_position
    for hv_pos in hv_positions:        
        # get files to process
        print('Get MOD13Q1 files to process')
        modis_file  = mm.select(product, year, day_of_year, hv_pos)
        
        if len(modis_file) == 0:
            print('File not exist')
            continue
        
        modis_file = modis_file[0]
        modis_files = mm.select(product, '*', day_of_year, hv_pos)
        
        # prepare dest file
        dest_file = pathlib.Path(modis_file).name
        dest_file = dest_file[:dest_file.rfind('.')] + '-vhi_drought.hdf'
        dest_file = pathlib.Path(drought_save_path) / dest_file
        
        if dest_file.exists():
            print('%s : result already exists, skip.'%hv_pos)
            continue
        
        dest_file = str(dest_file) 
        
        ndvi_min, ndvi_mean, ndvi_max, ndvi_stat_mask, ndvi_mask_count = get_ndvi_stat(modis_files)
        bt_min, bt_mean, bt_max, bt_stat_mask, bt_mask_count = get_bt_stat(modis_files, k1, k2)
        vhi_min, vhi_mean, vhi_max, vhi_stat_mask, vhi_mask_count = get_vhi_stat(modis_files, 
                                                                                 k1, k2, k_tci, k_vci, 
                                                                                 ndvi_min, ndvi_max, ndvi_stat_mask, 
                                                                                 bt_min, bt_max, bt_stat_mask)
        
        print('Calculate vhi')
        vhi, vhi_mask = get_vhi(modis_file, k1, k2, k_tci, k_vci, ndvi_min, ndvi_max, ndvi_stat_mask, bt_min, bt_max, bt_stat_mask)
                
        print('Calculate vhi drought and grade')
        vhi_drought, vhi_drought_mask = calc_vhi_drought(vhi, vhi_mask, vhi_mean, vhi_stat_mask)
        vhi_drought_grade, vhi_drought_grade_mask = calc_vhi_drought_grade(vhi, vhi_mask, grades, vhi_min, vhi_mean, vhi_stat_mask)

        if filter_by_vi_quality:
            hdf_obj = pyhdf.SD.SD(modis_file)
            ndvi_quality, ndvi_quality_mask = get_dataframe(hdf_obj, '250m 16 days VI Quality')
            ndvi_quality = (ndvi_quality & 0x03) == 0 # NDVI quality
            
            vhi_drought            *= ndvi_quality
            vhi_drought_mask       *= ndvi_quality
            vhi_drought_grade      *= ndvi_quality
            vhi_drought_grade_mask *= ndvi_quality
            hdf_obj.end()
        
        # save result
        print('Saving drought value and grade ...')
        
        vhi_drought *= 1000
        HdfMan.replace_hdf_dataframe(modis_file, dest_file,
                                     [vhi_drought, vhi_drought_grade],
                                     [vhi_drought_mask, vhi_drought_grade_mask],
                                     ['250m 16 days NDVI', '250m 16 days EVI'], 
                                     [-2000, -2000])
        
        print('---- done. ----\n')
    return


'''
if sys.platform == 'win32':
    modis_path = r'D:\CUIT-NewFarman\Data\MOD13Q1\h27v05'
else:
    modis_path = '/data/yule/disk2/MOD13Q1'
modis_hdf_file_list_file = 'file_list.txt'
ListFiles.list_files_recursive_to_file(modis_path, '*.hdf', modis_hdf_file_list_file)

# init
print('Init input varables')
modis_hdf_file_list_file = 'file_list.txt'
product = 'MOD13Q1'
year = 2015
day_of_year = 177
hv_pos = 'h27v05'
k1 = 60.776
k2 = 1560.56
k_vci = 1
k_tci = 1
grades = 5

if sys.platform == 'win32':
    drought_save_path = r'vhi_drought_save_path'
else:
    drought_save_path = '/data/yule/disk2/ProgramData/vhi_drought_test'

filter_by_vi_quality = False

calculate_vhi_drought_with_coor(modis_hdf_file_list_file, 
                                 hv_pos, 
                                 year, 
                                 day_of_year,
                                 k1, k2, k_vci, k_tci,
                                 grades,                                 
                                 drought_save_path,
                                 filter_by_vi_quality)
#'''

#-----------------------------------------------------------------------------

if __name__ == '__main__':
    if len(sys.argv) != 12:
        print('\nUsage')
        print('-' * 40)
        print('    python VHIDrought.py modis_file_path hv_pos year day_of_year k1 k2 k_vci k_tci grades drought_save_path filter_by_vi_quality')
        print()
        print('ndvi drought value is saved to "250m 16 days NDVI"       layer of dest hdf file.')
        print('ndvi drought grade is saved to "250m 16 days VI Quality" layer of dest hdf file.')
        print()
    else: # calculate ndvi drought value        
        modis_hdf_file_list_file, hv_pos, year, day_of_year, k1, k2, k_vci, k_tci, grades, drought_path, filter_by_vi_quality = sys.argv[1:]        
        calculate_vhi_drought_with_coor(modis_hdf_file_list_file, hv_pos, year, day_of_year, k1, k2, k_vci, k_tci, grades, drought_path, filter_by_vi_quality)
